#ifndef LUA_CPP_STD_VECTOR_UI
#define LUA_CPP_STD_VECTOR_UI
#include <luaCPPSTD/config.h>
#include <vector>
#include <numeric>
#include <complex>
#include <algorithm>
#include <iostream>
#include <boost/bind.hpp>
//after run tolua, disable tolua_function
//tolua_function(tolua_S,".add",tolua_luaVector_std_vector_std__vector____add01);

/*//////

// tolua_begin
namespace std{ 
	class vector
	{
		TOLUA_TEMPLATE_BIND(TT,double,std::complex<double>,std::vector<double>)

		vector();
		vector(int);
		vector(int,TT);
			
		//const TT& operator[](int index) const;
        //TT& operator[](int index);

		tolua_outside TT vector_geti @ .geti(int idx);
		tolua_outside void vector_seti @ .seti(int idx,const TT& value);


		void push_back(const TT& _Val);
		void reserve(int);
		void clear();
        int size() const;	
		void resize(int);
		void resize(int,const TT&);
		bool empty() const;

		tolua_outside std::vector<TT> add @ .add( double );
		tolua_outside std::vector<TT> sub @ .sub( double );
		tolua_outside std::vector<TT> mul @ .mul( double );
		tolua_outside std::vector<TT> div @ .div( double );

		tolua_outside std::vector< std::complex< double > > add @ .add( const std::complex<double>& );
		tolua_outside std::vector< std::complex< double > > sub @ .sub( const std::complex<double>& );
		tolua_outside std::vector< std::complex< double > > mul @ .mul( const std::complex<double>& );
		tolua_outside std::vector< std::complex< double > > div @ .div( const std::complex<double>& );
	};


	double min( const vector<double>& );
	double max( const vector<double>& );
	double sum( const vector<double>& );
	void min_and_index( const vector<double>& vec, double* val, int* idx );
	void min_and_index_with_range( const vector<double>& vec, int min, int max, double* val, int* idx );
	void max_and_index( const vector<double>& vec, double* val, int* idx );
	void max_and_index_with_range( const vector<double>& vec, int min, int max, double* val, int* idx );

	double interpolation(const vector<double>& xvec, const vector<double>& yvec,double x);



	vector<double> log10( const vector<double>& );
	vector<std::complex<double> > log10( const vector<std::complex<double> >& );
	vector<double> log( const vector<double>& );
	vector<std::complex<double> > log( const vector<std::complex<double> >& );
	vector<double> pow( const vector<double>&, double );
	vector<std::complex<double> > pow( const vector<std::complex<double> >& ,double);
	vector<double> abs( const vector<double>& );
	vector<double> abs( const vector<std::complex<double> >& );
	vector<double> arg( const vector<std::complex<double> >& );
	vector<double> phase( const vector<std::complex<double> >& );
	vector<double> imag( const vector<std::complex<double> >& );
	vector<double> real( const vector<std::complex<double> >& );

	vector<double> sin( const vector<double>& );
	vector<double> cos( const vector<double>& );
	vector<double> sinh( const vector<double>& );
	vector<double> cosh( const vector<double>& );
	vector<double> asin( const vector<double>& );
	vector<double> acos( const vector<double>& );


	//std::complex<double> min( const vector<std::complex<double> >& );
	//std::complex<double> max( const vector<std::complex<double> >& );

	double abs( const std::complex<double>& );
	double arg( const std::complex<double>& );
	double real( const std::complex<double>& );
	double imag( const std::complex<double>& );
	std::complex<double> polar(double);
	std::complex<double> polar(double,double);

	std::complex<double> sqrt(std::complex<double>);
	std::complex<double> pow(std::complex<double>,double);


	std::complex<double> acos(std::complex<double>);
	std::complex<double> acosh(std::complex<double>);
	std::complex<double> asin(std::complex<double>);
	std::complex<double> asinh(std::complex<double>);
	std::complex<double> atan(std::complex<double>);
	std::complex<double> atanh(std::complex<double>);

	std::complex<double> cos(std::complex<double>);
	std::complex<double> cosh(std::complex<double>);
	std::complex<double> sin(std::complex<double>);
	std::complex<double> sinh(std::complex<double>);
	std::complex<double> tan(std::complex<double>);
	std::complex<double> tanh(std::complex<double>);

	std::complex<double> exp(std::complex<double>);
};

$[

-- Create a vector constructor
function NewVector(self)
    local cobj = std.vector:new()
    for k,v in pairs(self) do
    	if type(v)=="number" then
    	    cobj:push_back(v)
    	end
    end
    tolua.takeownership(cobj)
    return cobj
end

function NewComplexVector(self)
    local cobj = std.vector_std__complex_double___:new()
    for k,v in pairs(self) do
    	if type(v)=="number" then
    	    cobj:push_back(v)
    	end
    end
    tolua.takeownership(cobj)
    return cobj
end

function LinearSpace(min,max,num, arraytype)
	assert( type(num) == "number", "num is not set" );
	assert( num >= 1, "num is wrong" );
	local cobj;
	if arraytype == nil or arraytype ~= "complex"  then
		assert( type(min) == "number", "min is not set" );
		assert( type(max) == "number", "max is not set" );
		cobj = std.vector:new(num);
    end

    if arraytype == "complex" then
		cobj = std.vector_std__complex_double___(num)
		min = std.complex:new(min)
		max = std.complex:new(max)
    end

    if num == 1 then
		cobj[0] = min;
	else
		local delta = (max-min)/(num-1);
		for ii=0, cobj:size()-1 do
			cobj[ii] = min + ii * delta ;
		end
	end
	--tolua.takeownership(cobj);
	return cobj;
end

$]

// tolua_end
//*//////

namespace std{

	template<class TT>
	TT vector_geti(std::vector<TT>* a, int idx)
	{
		if(!a || idx<0 || idx >= a->size() )
			throw std::out_of_range("out of range: "+std::to_string(idx));
		return a->operator [](idx);
	}

	template<class TT>
	void vector_seti(std::vector<TT>* a, int idx,const TT& v)
	{
		if(!a || idx<0 || idx >= a->size() )
			throw std::out_of_range("out of range: "+std::to_string(idx));
		a->operator [](idx) = v;
	}

	template<class T> 
	T min( const vector<T>& vec )
	{
		return *std::min_element(vec.begin(),vec.end());
	}

	template<class T>
	void min_and_index( const vector<T>& vec, T& val, int& idx )
	{
		auto ele = std::min_element(vec.begin(),vec.end());
		val = *ele;
		idx = std::distance(vec.begin(),ele);
	}

	template<class T>
	void min_and_index_with_range( const vector<T>& vec, int min, int max, T& val, int& idx )
	{
		if( min > max ) std::swap(min,max);
		min = std::min( std::max( 0, min), (int)vec.size() );
		max = std::max( std::min( (int)vec.size(), max ), 0 );
		auto it1 = vec.begin()+min;
		auto it2 = vec.begin()+max;
		auto ele = std::min_element(it1,it2);
		val = *ele;
		idx = std::distance(vec.begin(),ele);
	}

	template<class T>
	void max_and_index( const vector<T>& vec, T& val, int& idx )
	{
		auto ele = std::max_element(vec.begin(),vec.end());
		val = *ele;
		idx = std::distance(vec.begin(),ele);
	}

	template<class T>
	void max_and_index_with_range( const vector<T>& vec, int min, int max, T& val, int& idx )
	{
		if( min > max ) std::swap(min,max);
		min = std::min( std::max( 0, min), (int)vec.size() );
		max = std::max( std::min( (int)vec.size(), max ), 0 );
		auto it1 = vec.begin()+min;
		auto it2 = vec.begin()+max;
		auto ele = std::max_element(it1,it2);
		val = *ele;
		idx = std::distance(vec.begin(),ele);
	}

	template<class T>
	T interpolation(const std::vector<T>& xarr, const std::vector<T>& yarr,T x)
	{
		auto pos = std::find_if(xarr.begin(),xarr.end(),boost::bind(std::greater<T>(),_1,x));
		int npos = std::distance(xarr.begin(),pos);
		int nlow = std::max(npos-1,0);
		int nhigh = std::min((int)xarr.size()-1,npos);
		if( npos == 0) return yarr[0];
		if( npos+1 == xarr.size() ) return yarr[xarr.size()-1];
		T res = yarr[npos-1] + (x-xarr[npos-1])/(xarr[npos]-xarr[npos-1])*(yarr[npos]-yarr[npos-1]);
		return res;
	}

	template<class T> 
	T max( const vector<T>& vec )
	{
		return *std::max_element(vec.begin(),vec.end());		
	}
	template<class T>
	T sum( const vector<T>& vec )
	{
		T res(0);
		return std::accumulate(vec.begin(),vec.end(),res);
	}

	template<class T>
	vector<T> log10( const vector<T>& vec )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::log10(res[ii]);
		return res;
	}

	template<class T>
	vector<T> log( const vector<T>& vec )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::log(res[ii]);
		return res;
	}

	template<class T>
	vector<T> sin( const vector<T>& vec )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::sin(res[ii]);
		return res;
	}

	template<class T>
	vector<T> cos( const vector<T>& vec )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::cos(res[ii]);
		return res;
	}

	template<class T>
	vector<T> asin( const vector<T>& vec )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::asin(res[ii]);
		return res;
	}

	template<class T>
	vector<T> acos( const vector<T>& vec )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::acos(res[ii]);
		return res;
	}

	template<class T>
	vector<T> sinh( const vector<T>& vec )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::sinh(res[ii]);
		return res;
	}

	template<class T>
	vector<T> cosh( const vector<T>& vec )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::cosh(res[ii]);
		return res;
	}



	template<class T>
	vector<T> pow( const vector<T>& vec, double e )
	{
		vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::pow(res[ii],e);
		return res;
	}

	template<class T>
	vector<double> abs( const vector<T>& vec )
	{
		vector<double> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::abs(vec[ii]);
		return res;
	}

	template<class T>
	vector<double> arg( const vector<T>& vec )
	{
		vector<double> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::arg(vec[ii]);
		return res;
	}
	template<class T>
	vector<double> phase( const vector<T>& vec )
	{
		vector<double> res(vec.size());
		double phs0;
		for( int ii=0; ii < res.size(); ii++)
		{
			if( ii == 0 )
			{
				phs0 = std::arg( vec[ii] );
			}
			else
			{
				phs0 = phs0 + std::arg(vec[ii]/vec[ii-1]);
			}
			res[ii] = phs0;
		}
		return res;
	}

	template<class T>
	vector<double> real( const vector<T>& vec )
	{
		vector<double> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::real(vec[ii]);
		return res;
	}

	template<class T>
	vector<double> imag( const vector<T>& vec )
	{
		vector<double> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::imag(vec[ii]);
		return res;
	}


	template<class TT>
		std::vector<TT> add(std::vector<TT>* a, double b);
	template<class TT>
		std::vector<TT> sub(std::vector<TT>* a, double b);
	template<class TT>
		std::vector<TT> mul(std::vector<TT>* a, double b);
	template<class TT>
		std::vector<TT> div(std::vector<TT>* a, double b);

	std::vector<int> add(std::vector<int>* a, int b);
	std::vector<int> sub(std::vector<int>* a, int b);
	std::vector<int> mul(std::vector<int>* a, int b);
	std::vector<int> div(std::vector<int>* a, int b);

	template<class TT>
	std::vector< std::complex<double> > add(std::vector<TT>* a,  const std::complex<double>& b )
	{
		std::vector< std::complex<double> > res( a->size() );
		for( int ii=0; ii < res.size(); ii ++ )
		{
			res[ii] = (*a)[ii] + b;
		}
		return res;
	}

	template<class TT>
	std::vector< std::complex<double> > sub(std::vector<TT>* a,  const std::complex<double>& b )
	{
		std::vector< std::complex<double> > res( a->size() );
		for( int ii=0; ii < res.size(); ii ++ )
		{
			res[ii] = (*a)[ii] - b;
		}
		return res;
	}

	template<class TT>
	std::vector< std::complex<double> > mul(std::vector<TT>* a,  const std::complex<double>& b )
	{
		std::vector< std::complex<double> > res( a->size() );
		for( int ii=0; ii < res.size(); ii ++ )
		{
			res[ii] = (*a)[ii] * b;
		}
		return res;
	}

	template<class TT>
	std::vector< std::complex<double> > div(std::vector<TT>* a,  const std::complex<double>& b )
	{
		std::vector< std::complex<double> > res( a->size() );
		for( int ii=0; ii < res.size(); ii ++ )
		{
			res[ii] = (*a)[ii] / b;
		}
		return res;
	}
};

#endif
